Model of random packings of different size balls 
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We develop a model to describe the properties of random assemblies of polydisperse hard spheres. 
We show that the key features to describe the system are (i) the dependence between the free volume 
of a sphere and the various coordination numbers between the species, and (ii) the dependence of the 
coordination numbers with the concentration of species; quantities that are calculated analytically. 
The model predicts the density of random close packing and random loose packing of polydisperse 
systems for a given distribution of ball size and describes packings for any interparticle friction 
coefficient. The formalism allows to determine the optimal packing over different distributions and 
may help to treat packing problems of non-spherical particles which are notoriously difficult to solve. 



Understanding the basic properties of spheres pack- 
ings is a major challenge since this problem may provide 
valuable knowledge regarding low temperature phases in 
condensed matter physics The canonical example is 
perhaps the monodisperse sphere packing problem. It 
has been mathematically proven that the optimum way 
to arrange monodisperse spheres is the face-centered cu- 
bic lattice; a problem that has been solved recently by 
Hales, ~ 400 years after the famous Kepler conjecture 
on the issue. On the other hand, it is commonly ob- 
served that packings arrange in a random fashion at a 
lower density state called random close packing or RCP 
Q . Furthermore, packings are mechanically stable up to 
an even lower limit called random loose packing, RLP. 

In parallel with the large literature dealing with 
monodisperse sphere packings, a large body of experi- 
mental, theoretical and numerical work has been devoted 
to the analysis of polydisperse systems; the interest aris- 
ing due to the simple fact that polydispersivity is om- 
nipresent in most realistic systems and industrial appli- 
cations While previous approaches have focused on 
frictionless packings, an integrated analytical approach 
that brings together different observations for all pack- 
ings from RLP to RCP and for any friction or coordi- 
nation number is still lacking. Based on our previous 
statistical mechanics approach Q, here we build such a 
framework. 

We show that the key aspect to solve this problem is 
the dependence of the various coordination numbers be- 
tween the different species and the concentration of the 
species. This is calculated here and shown to agree well 
with computer simulations. This result is then incorpo- 
rated into a statistical theory of volume fluctuations as 
in Q which calculates the free volume of a particle in 
terms of the coordination number. The main result is 
the prediction of the RLP and RCP limiting densities for 
a given distribution of ball sizes as well as the prediction 
of densities for any packing in between those limits. The 
formalism allows for a determination of the best pack- 
ing fraction in terms of different distribution of ball sizes 
with specified constraints, as we show with a simple ex- 
ample. We discuss possible generalization of the method 
to solve more difficult problems like the phase behavior 



of systems of non-spherical particles like rods or sphero- 
cylinders in any dimensions; problems of long-standing 
history in condensed matter @. 

Recent theoretical advances [|[ allow for the prediction 
of the density of RCP and RLP for equal-size ball pack- 
ings using a relation between the average volume and 
the geometrical coordination number. Following this ap- 
proach, we here describe long-range spatial correlations 
through a mean-field background term. This approxi- 
mation makes the problem amenable to analytic calcu- 
lations, and is shown to describe well our simulation re- 
sults. An explicit inclusion of such correlations is possible 
in our framework, but severely complicates any solution 
attempts. Thus, we believe that the present approach is 
accurate enough for many important properties, such as 
the volume fraction calculation. 

The above theoretical framework will guide the present 
formalism for polydisperse systems. We first treat the 
case of binary mixtures of hard spheres of radius R\ and 
R 2 > Ri in 3d and then generalize the problem to any 
distribution and dimension. 

Calculation of z^. — The key quantity to calculate is 
Zij , denoting the mean number of contacts of a ball of ra- 
dius Ri with a ball of radius Rj , versus the concentration 
of one of the species. We need a formula for Zij tlS el func- 
tion of (z,x, J^), the later being the size ratio, x is the 
fraction of small balls in the packing x = Nx/(N\ + N 2 ) 
with Ni the number of i balls and z is the global geomet- 
rical coordination number averaged over all the particles: 
z = xz\ + (1 — x)z2 where Zi is the average coordination 
of each species. 

The coordinations are determined by three equations: 

Zi = z a + z i2 , (i = 1,2), xz 12 = (1 - x)z 21 . (1) 

We assume that these coordinations are inversely pro- 
portional to the average solid angle extended by con- 
tacting balls Ri and R 2 . The average solid angles are 
denoted (S° cc ) and are calculated in terms of the solid 
angle that a ball Rj occupy on a ball Ri according to 



2 



qOCC 




FIG. 1: (a) Occupied surface and free surface, (b) Voronoi cell 
for polydisperse balls. Plots are in 2d for easier visualization. 
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FIG. 2: Comparison between theory and numerical simula- 
tions for (a) Zij, Eq. (g}, and (b) P>(c), Eq. (fTTT) . 




(S° cc ) = xSff c + (1 - x)Sf£ c with (see Fig. Q}t) 
S° J cc /2n z 

Thus, (S° cc ) represents the mean occupied surface on 
a i ball weighted by the concentrations x and (1 — x). 
This represents an approximation since the real weights 
are zn/zi and z^/^i, respectively. Then, Zi oc 1/(S° CC ) 
leading to the following normalizations: 



the classical Voronoi polyhedron defined by the center of 
the particle, one should consider all the points which are 
closer to the surface of a given particle. Such a construc- 
tion is called the Voronoi S region and tiles a system of 
nonspherical convex particles and polydisperse systems 
as can be seen in Fig. []J>. Following this approach we 
calculate the average volume of a polydisperse Voronoi 
cell, denoted W . The volume fraction is given by (j) = -j^-, 
where V g = ^-{xR\ + (1 — x)R 2 ) is the mean volume of a 
ball. We first find the analytical formula of the Voronoi S 
region. The boundary of the Voronoi cell in the direction 
s of a i ball next to a j ball at position ry is (Fig. [1}d): 



z + (l-z)|fp' 
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Thus, the system of Eqs. $Q is reduced to a system 

of three equations for four unknowns z ir To close the where s and hj are unitary. Thus, the boundary of a 

system we assume proportional laws and deduce from Voronoi cell of a ball i in the direction s is the minimum 

Zi by considering that z l} is proportional to the number f over a u the particles j for any ^ (s) > 0. This 



of contacts of the i balls times the number of contacts of 
the j balls: 



zn = A(z 1 x){z 1 x), z 12 
z 21 = Bz 2 (l - x)(z 1 x), z 22 



A( Zl x)z 2 (l-x), (3) 
Bz 2 (l - x)z 2 {\ - x). 



leads to 
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Using the first equation in ([lj we find the constants A 
and B, leading to the solution: 
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The volume of the cell of the ball i is then given by 
Wi = i § li (s) 3 ds. We define the orientational Voronoi 
volume, Wf, along the direction s by Wi = § Wfds = 
(W^s. This leads to 



(4) 



Figure [5^, compares this solution to numerical simula- 
tions of Hertz packings jammed at RCP [5j for Ri/R 2 = 
1.4 and z = 6. We find that the formulae are very accu- 
rate for size ratios below 1.5 and present small deviations 
up to size ratio 2. The results are also in agreement with 
II- 

The Voronoi cell.— The common way to divide a 
system into volumes associated with each particle is the 
Voronoi tessellation. The Voronoi cell for monodisperse 
particles [|| is composed by all the points nearest to the 
center of the ball than to any other ball. This definition 
has been extended in Q to the case of polydisperse sys- 
tems and non-spherical particles: instead of considering 



Wf = - 



Tij .S> 



rj j -{R i -R j Y 

T *- R j ^nMj.s) - (Ri — Rj) 



(7) 



This definition leads to the results of |5j when R\ = R 2 . 
Since the system is isotropic, the mean Voronoi volume 
can be calculated as: 



(8) 



Calculation of the mean Voronoi volume. — Hav- 
ing calculated the Voronoi cell exactly in Eq. ([7]), we now 
proceed to develop a probability theory of volume fluc- 
tuations in the spirit of the quasiparticle approximation 
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used in Q to obtain the mean Voronoi volume. For a 
given ball i, the calculation of W* reduces to finding the 
ball j* that minimizes la(s). We call j* the Voronoi 
ball for the ball i. 



Hj(s)- We call j 
We consider r = r,,* , cost/ = s ■ r. 



and c = 2L 



We have c 
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Therefore, 



"V ■ ,,c iiclv ^ ° rcosO-iRi-Rj) ■ 

we just need to compute the inverse cumulative distribu- 
tion function, denoted P>(c), to find all the balls j with 
Uj > I , and thus not contributing to the Voronoi volume 
of the ball i. The average Voronoi volume is then given 
by the expression 
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dc 2 
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c)dc. 



We calculate the mean Voronoi volume for the balls 
of radius Pi and P2 separately and then average them. 
We denote P 1 , (c) and P 2 , (c) the inverse cumulative dis- 



tributions respectively, and W 

*)U°°- 2d2 



P 2 (c)cfc, and therefore, 



P > (c) = xP*(c) + (l-x)I#(c). 



(10) 



Calculation of P>(c). — There are three salient steps 
in the calculation of P>(c): (i) The separation of P>(c) 
following Eq. (fTU)) . (m) The separation of each term 
P>(c), i = 1,2, into two contributions: a term taking into 
account the contact spheres, P> c '(c), and a bulk term, 
P> (c). The contact term clearly depends on Zij. The 
bulk term averages over all spatial correlations of non- 
contact particles and, without significant loss of accuracy 
as shown below, we assume that it only depends on the 
average value of W. In principle, it is possible to use a 
more realistic form for this term, but this would render 
the problem practically unsolvable. (Hi) The separation 
of Pf (c) and Pi? (c) into two terms P> jC (c) and P> jS (c) , 
j = 1,2, for each species. 

An assumption of the theory (to be tested a posteriori 
with computer simulations) is that all of these terms are 
independent. Thus 



P>(c) 



xPl lc (c)Pl 2C (c)Pl 1B (c)Pl 2B (c) + 
(1 - x)Pl lc {c)Pl 2C {c)Pl 1B {c)Pl 2B {c). 



(11) 



Also, we work in the limit of large number of parti- 
cles leading to Boltzmann-like exponential distributions 



for each P 
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and Pi? Q. Four important quantities 
are then defined: (i) V*j{c) and (ii) <Sy(c): the excluded 
volume and surface on the ball, respectively, where no 
center of a ball j can be located for a given ball i and for 
a given c. (in) pf. the mean number of balls j by the left 
free volume, (iv) p\y the mean number of balls j by the 
left free surface on a ball i. These considerations lead to: 



P» B (c) 



exp 



(-PjV*{c) 
P!f c (c)=exp(-^5*(c) 



(12) 



Next, we calculate these four quantities. To simplify 



we denote / = Ri + Rj, k 
function. We obtain: 



Ri — Rj and the step- 
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(9) where x\ = x and x% = (1 — x). The fourth quantity, 



pfj, is the most difficult to calculate. In terms of the oc- 
cupied areas Eqs. © we have pfj = 4T _ 2 . 1 g. Z c 'c_ z . 2 s°cc ■ 
However, for the contact terms, the analogy with the 
Boltzmann-like exponential distribution of volumes is far 
from being exact. This is because this form is based 
on the large number limit which in the case of contact- 
ing balls is no more than around 6. Therefore, the ex- 
ponential distribution is used as an Ansatz with p'- a 
variational parameter as in @. We denote (5|f e ) the 
mean solid angle of the gaps left between the j contact- 
ing neighbors of a i ball (Fig. [Th). We obtain: 



(^ ee ) 
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d(l-P>(c)) 
dc 



dc 



(14) 
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Then, we perform numerical simulations to find 
{S{j ee ). We randomly generate balls of radius Ri and Rj 
with the proportion zn/zi and z i2 / z i respectively around 
a ball of radius Ri and then evaluate the mean free sur- 
face (S-j £e ) and its inverse to obtain p\-. We find 
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(15) 
■) 2 )), 



which is a generalization of the results of @ to polydis- 
perse systems. 

Using Eqs. (HI]), (H3J) and ^ into ([II]), P>(c) can be 
calculated by solving the equations numerically. Figure 
^Bp depicts the comparison of the theoretical results of the 
probability of Voronoi volumes P>(c) with computer gen- 
erated Hertzian packings with z = 6 for x = 0.5 at RCP. 
The calculated distribution is quite accurate for most of 
the range except for small deviations at high values of 
c, which however, do not affect much the value of the 
average (c 3 ) in Eq. ©. This shows that our mean-field 
approximation already captures the main contribution of 
the probability distribution function P>(c). 
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Calculation of W. — The above considerations lead 
to a final form to calculate W using Eq. (fTTT) into ©: 

/"°° / \ 

(16) 

Notice that Pj-(W) depends on JV, Eq. (|L3"1) . and pfj(zij) 
depends on the Zij, Eq. (TT51) . which in turn depends on 
the concentration a; and z through Eq. ([4]). Therefore 
Eq. (ITB1) is a self-consistent equation to obtain W(^, x), 
for a given Ri/B.2. A numerical integration of Eq. (fT6|) 
is performed versus a; for a given value of z. By consid- 
ering the isostatic limits of z = 6 and z = 4 we predict 
the limits of RCP and RLP at zero friction and infinite 
friction between the spheres, respectively @. The results 
for the volume fraction at RCP versus x are depicted in 
Fig. [SJi which also compares the results to numerically 
generated packings of Hertz spheres [f|. We see a very 
good agreement indicating that the theory captures well 
the behavior of polydisperse packings and that the ap- 
proximations used are reasonable. For size ratios larger 
than 2 deviations are found indicating the limit of valid- 
ity of the approach. For any other value of interparticle 
friction between and oo, the density is obtained by set- 
ting z between the limiting isostatic values of 6 and 4, 
respectively. The resulting volume fraction is shown in 
Fig. Eb- The result for RLP for a given x is a new pre- 
diction as this problem has not been investigated before. 
Our results promote new experiments to test the RLP 
limit of polydisperse systems shown in Fig. [3Jd. 

The formalism can be extended to consider any distri- 
bution of sphere radius. The main modification is that 
all the sums over the radius are replaced by integrations 
over the desired distribution of radius P(R) (the binary 
case corresponds to two delta- functions at R\ and i?2)- 

And all the quantities are calculated for balls of inter- 
nal radius r and external R and x and (1 — x) are replaced 
by P(r)dr and P(R)dR, respectively, including the coor- 
dinations (see Supplementary Information B). This re- 
sult allows to explore the space of radius distributions in 
search of the optimal packings for a given polydispersiv- 
ity. This analysis could be of industrial interest if one 
wishes to optimize the packing fraction by introducing 
different species in the mixture. 

We calculate the volume fraction for several distribu- 
tions P{R) constraint by ball radius in the range [1,2], in 
search of the optimal packing. We calculate the volume 
fraction for various P(R) ranging from uniform to two- 
peaked Gaussian distributions of varied widths. We find 
that the more small balls we add the better the packing 
until a certain point where the volume fraction starts to 
decrease. This maximum can be rationalized assuming 
that the small balls always fill the gaps between the large 
ones as long as there are enough large balls. Further ex- 
tensions of the theory to any dimension can be performed 
by replacing 3 by d in Eq. (JT]) and developing a theory 
of volume fluctuation in d-dimensions. We notice that 



many of the appro ximations employed in 3 d may become 
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FIG. 3: (a) Comparison between theory and numerical sim- 
ulations of Hertzian packings at RCP, i.e. z — 6 versus x 
for different values of R1/R2 as indicated. Error bars are std 
over 10 realizations of the packings with 10,000 balls, (b) 
Three dimensional surface plot of cj> as a function of z and x 
for R1/R2 = 1.5. The numerical results at RCP and RLP are 
provided in Supplementary Information A. 



exact for large d, thus we expect improved results in the 
mean field limit of infinite dimensions. 

The method allows to treat more difficult problems. 
For instance, the prediction of the volume fraction of 
a jammed system of non-spherical particles is a long- 
standing problem. Theoretical predictions of Onsager 
@ are valid for large aspect ratios, like elongated rods. 
Experiments however, find interesting new physics for 
small aspect ratios. In this respect, the present polydis- 
perse theory could be mapped to the problem of ellip- 
soids, spherocylinders or rods. A Voronoi cell needs to 
be calculated as a function of the angles defining the ori- 
entation of the non-spherical particles in analogy of the 
calculation between two particles of different radii. The 
integration over P(r)dr in Eq. (|16[) is then replaced by in- 
tegration over weighted oricntational angles. The above 
analysis can also be extended to dimensions beyond three 
8] . Although many of the appproximations should work 
better in higher dimensions, some of the hypotheses (for 
example, the contact term ansatz) need to be reassessed. 
Thus, higher dimensions studies cannot be addressed as 
trivial extensions and need to be handled with care. 

In summary, a theoretical framework is presented that 
predicts the RLP and RCP limits of a system of polydis- 
perse spheres and brings together distinct results into a 
common theoretical framework. The formalism has the 
potential to solve other problems in condensed matter 
physics such as the mixing and phase behavior of sys- 
tems of hard particles of different shapes and size. 
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Appendix A: Table of volume fraction values predicted by the theory. 



x 





0.1 


0.2 


0.3 


0.4 


0.5 


0.6 


0.7 


0.8 


0.9 


1.3 


0.5359 


0.5376 


0.5392 


0.5405 


0.5416 


0.5423 


0.5426 


0.5423 


0.5413 


0.5393 


1.5 


0.5359 


0.5393 


0.5426 


0.5456 


0.5482 


0.5503 


0.5516 


0.5517 


0.5499 


0.5453 


1.7 


0.5359 


0.5409 


0.5457 


0.5506 


0.5550 


0.5588 


0.5616 


0.5628 


0.5611 


0.5359 



TABLE I: Volume fraction for z = 4, RLP. 








0.1 


0.2 


0.3 


0.4 


0.5 


0.6 


0.7 


0.8 


0.9 


1.3 


0.6340 


0.6356 


0.6370 


0.6382 


0.6392 


0.6399 


0.6402 


0.6399 


0.6389 


0.6340 


1.5 


0.6340 


0.6372 


0.6401 


0.6423 


0.6453 


0.6472 


0.6484 


0.6484 


0.6468 


0.6426 


1.7 


0.6340 


0.6386 


0.6431 


0.6475 


0.6515 


0.6549 


0.6575 


0.6586 


0.6571 


0.6506 



TABLE II: Volume fraction for 2 = 6, RCP. 



Appendix B: Distribution of sphere radius. 



According to the theory, the average Voronoi volume for a packing with a distribution of radius P(r), is given by 
the following self-consistent equation : 



J drP{r)Jc 2 cxp ( J dR(-p s rR S* rR (c) - p R {W)V r * R {c))) dc, (Bl) 



where the different quantities are calculated as follow: 

VMc) = J 8(c - r ^/_\ )A- 3 = 2*(-i«c - kf - C) + l((c - k)° - l') + (g - |)((c - kf - I 2 )), 

To simplify we denoted I = r + R, k = r — R and <d the step-function. 
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